今日やること!
ggplotでベクター形式のGISデータを可視化する
ggplotでラスター形式のGISデータを可視化する
2016年のアメリカ大統領選挙の、州ごとの投票率および各得票数に関する変数のデータを可視化する。
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(socviz)
election %>%
select(fips, state, total_vote, r_points, pct_trump, party, census) %>%
slice_sample(n = 5)
## # A tibble: 5 × 7
## fips state total_vote r_points pct_trump party census
## <dbl> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 17 Illinois 5594825 -16.9 38.4 Democratic Midwest
## 2 11 District of Columbia 311268 -86.8 4.09 Democratic South
## 3 42 Pennsylvania 6166729 0.710 48.2 Republican Northeast
## 4 26 Michigan 4824542 0.220 47.2 Republican Midwest
## 5 48 Texas 8993166 8.98 52.1 Republican South
FIPS (federal information processin standard) コードは、州・準州に割り当てられる連邦情報処理標準による番号。
party_colors <- c('#2E74C0', '#CB454A')
p0 <- ggplot(data = subset(election, st %nin% 'DC'),
mapping = aes(x = r_points,
y = reorder(state, r_points),
color = party)) +
geom_vline(xintercept = 0, color = 'grey30') +
geom_point(size = 2)
p0
p1 <- p0 +
scale_color_manual(values = party_colors) +
scale_x_continuous(breaks = c(-30, -20, -10, 0, 10, 20, 30, 40),
labels = c('30\n(Clinton)', '20', '10', '0',
'10', '20', '30', '40\n(Trump)')) +
facet_wrap(~census, ncol = 1, scales = 'free_y') +
guides(color = 'none') +
labs(x = 'Point Margin', y = '') +
theme(axis.text = element_text(size = 8)) +
theme_bw()
p1
「空間データを必ずしも空間構造として表現しなくても良い」
library(maps)
##
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
##
## map
us_states <- map_data('state')
class(us_states)
## [1] "data.frame"
str(us_states)
## 'data.frame': 15537 obs. of 6 variables:
## $ long : num -87.5 -87.5 -87.5 -87.5 -87.6 ...
## $ lat : num 30.4 30.4 30.4 30.3 30.3 ...
## $ group : num 1 1 1 1 1 1 1 1 1 1 ...
## $ order : int 1 2 3 4 5 6 7 8 9 10 ...
## $ region : chr "alabama" "alabama" "alabama" "alabama" ...
## $ subregion: chr NA NA NA NA ...
head(us_states)
## long lat group order region subregion
## 1 -87.46201 30.38968 1 1 alabama <NA>
## 2 -87.48493 30.37249 1 2 alabama <NA>
## 3 -87.52503 30.37249 1 3 alabama <NA>
## 4 -87.53076 30.33239 1 4 alabama <NA>
## 5 -87.57087 30.32665 1 5 alabama <NA>
## 6 -87.58806 30.32665 1 6 alabama <NA>
dim(us_states)
## [1] 15537 6
まずは白地図を書く。
p <- ggplot(data = us_states,
mapping = aes(x = long, y = lat,
group = group)) +
geom_polygon(fill = 'white', color = 'black')
p
州ごとに色分けしたコロプレス図を描いてみる。
p <- ggplot(data = us_states,
mapping = aes(x = long, y = lat,
group = group,
fill = region)) +
geom_polygon(color = 'grey90', size = 0.1) +
guides(fill = 'none')
p
デフォルトではメルカトル図法で書かれているが、広域のGISデータの描画では歪みが生じる。アルベルス正積円錐図法を使うことで、面積の比率が正しくなる。
p <- ggplot(data = us_states,
mapping = aes(x = long, y = lat,
group = group,
fill = region)) +
geom_polygon(color = 'grey90', size = 0.1) +
coord_map(projection = 'albers', lat0 = 39, lat1 = 45) + # lat0, lat1はアメリカの場合のデフォルト
guides(fill = 'none')
p
先程の州ごとのデータに、州名をキーにして、選挙関連の統計情報をdplyr::left_join()関数で結合する。
head(election)
## # A tibble: 6 × 22
## state st fips total_vote vote_margin winner party pct_margin r_points
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl>
## 1 Alabama AL 1 2123372 588708 Trump Repu… 0.277 27.7
## 2 Alaska AK 2 318608 46933 Trump Repu… 0.147 14.7
## 3 Arizona AZ 4 2604657 91234 Trump Repu… 0.035 3.5
## 4 Arkansas AR 5 1130635 304378 Trump Repu… 0.269 26.9
## 5 California CA 6 14237893 4269978 Clint… Demo… 0.300 -30.0
## 6 Colorado CO 8 2780247 136386 Clint… Demo… 0.0491 -4.91
## # … with 13 more variables: d_points <dbl>, pct_clinton <dbl>, pct_trump <dbl>,
## # pct_johnson <dbl>, pct_other <dbl>, clinton_vote <dbl>, trump_vote <dbl>,
## # johnson_vote <dbl>, other_vote <dbl>, ev_dem <dbl>, ev_rep <dbl>,
## # ev_oth <dbl>, census <chr>
head(us_states)
## long lat group order region subregion
## 1 -87.46201 30.38968 1 1 alabama <NA>
## 2 -87.48493 30.37249 1 2 alabama <NA>
## 3 -87.52503 30.37249 1 3 alabama <NA>
## 4 -87.53076 30.33239 1 4 alabama <NA>
## 5 -87.57087 30.32665 1 5 alabama <NA>
## 6 -87.58806 30.32665 1 6 alabama <NA>
# electionの州名を、us_statesと同様に全て小文字に変換する
# <-- よくある変換。そしてエラーが出やすい
election$region <- tolower(election$state)
# region列をキーにして、us_statesデータにelectionデータを結合する。
# left_join(): 左のus_statesデータを優先して、electionデータを結合する
# right_join(): 右のelectionデータを優先して、us_statesデータを結合する
# inner_join(): us_statesとelectionデータに共通する列だけ結合する
# full_join(): どちらか一方に含まれる列が結合される。存在しない列はNAになる
us_states_elec <- left_join(us_states, election, by = 'region')
head(us_states_elec)
## long lat group order region subregion state st fips total_vote
## 1 -87.46201 30.38968 1 1 alabama <NA> Alabama AL 1 2123372
## 2 -87.48493 30.37249 1 2 alabama <NA> Alabama AL 1 2123372
## 3 -87.52503 30.37249 1 3 alabama <NA> Alabama AL 1 2123372
## 4 -87.53076 30.33239 1 4 alabama <NA> Alabama AL 1 2123372
## 5 -87.57087 30.32665 1 5 alabama <NA> Alabama AL 1 2123372
## 6 -87.58806 30.32665 1 6 alabama <NA> Alabama AL 1 2123372
## vote_margin winner party pct_margin r_points d_points pct_clinton
## 1 588708 Trump Republican 0.2773 27.72 -27.72 34.36
## 2 588708 Trump Republican 0.2773 27.72 -27.72 34.36
## 3 588708 Trump Republican 0.2773 27.72 -27.72 34.36
## 4 588708 Trump Republican 0.2773 27.72 -27.72 34.36
## 5 588708 Trump Republican 0.2773 27.72 -27.72 34.36
## 6 588708 Trump Republican 0.2773 27.72 -27.72 34.36
## pct_trump pct_johnson pct_other clinton_vote trump_vote johnson_vote
## 1 62.08 2.09 1.46 729547 1318255 44467
## 2 62.08 2.09 1.46 729547 1318255 44467
## 3 62.08 2.09 1.46 729547 1318255 44467
## 4 62.08 2.09 1.46 729547 1318255 44467
## 5 62.08 2.09 1.46 729547 1318255 44467
## 6 62.08 2.09 1.46 729547 1318255 44467
## other_vote ev_dem ev_rep ev_oth census
## 1 31103 0 9 0 South
## 2 31103 0 9 0 South
## 3 31103 0 9 0 South
## 4 31103 0 9 0 South
## 5 31103 0 9 0 South
## 6 31103 0 9 0 South
選挙結果で共産党と民主党のどちらが勝利したのかを可視化してみる。
p <- ggplot(data = us_states_elec,
aes(x = long, y = lat,
group = group,
fill = party)) +
geom_polygon(color = 'grey90', size = 0.1) +
coord_map(projection = 'albers', lat0 = 39, lat1 = 45)
p
最初の作図の色と揃えてみる。ついでに見た目をcowplot::theme_map()で調整
library(cowplot)
p <- ggplot(data = us_states_elec,
aes(x = long, y = lat,
group = group,
fill = party)) +
geom_polygon(color = 'grey90', size = 0.1) +
coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
scale_fill_manual(values = party_colors) +
labs(title = 'Election Results 2016',
fill = NULL) +
theme_map() +
theme(legend.position = 'bottom')
p
ドナルド・トランプ氏の得票率を示してみる
p <- ggplot(data = us_states_elec,
aes(x = long, y = lat,
group = group,
fill = pct_trump)) +
geom_polygon(color = 'grey90', size = 0.1) +
coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
scale_fill_gradient(low = 'white',
high = '#CB454A') +
labs(title = 'Trump vote',
fill = 'Percent',
fill = NULL) +
theme_map() +
theme(legend.position = 'bottom')
p
対立する2党への投票の結果を、中間点から分岐する配色で示してみる。
パターン1. 中間色を白にする。
p0 <- ggplot(data = us_states_elec,
aes(x = long, y = lat,
group = group,
fill = d_points)) +
geom_polygon(color = 'grey90', size = 0.1) +
coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
labs(title = 'Winning margins',
fill = 'Percent') +
theme_map() +
theme(legend.position = 'bottom')
p1 <- p0 +
scale_fill_gradient2()
p1
パターン2: 中間を紫
p2 <- p0 +
scale_fill_gradient2(low = party_colors[2],
mid = scales::muted('purple'),
high = party_colors[1],
breaks = c(-25, 0, 25, 50, 75))
p2
パターン3: 民主党が強い基盤を持っているワシントンD.C.を除外して可視化
p3 <- ggplot(data = filter(us_states_elec, region %nin% 'district of columbia'),
aes(x = long, y = lat,
group = group,
fill = d_points)) +
geom_polygon(color = 'grey90', size = 0.1) +
coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
scale_fill_gradient2(low = party_colors[2],
mid = scales::muted('purple'),
high = party_colors[1],
breaks = c(-25, 0, 25, 50, 75)) +
labs(title = 'Winning margins',
fill = 'Percent') +
theme_map() +
theme(legend.position = 'bottom')
p3
コロプレス図は、マッピングしている変数を部分的にしか表現できない。地図+地図に頼らない可視化はセットで作ると良い。
ポリゴンデータ
ラスタデータ
ポリゴンデータ
ラスタデータ
terraパッケージ
starsパッケージ
https://notchained.hatenablog.com/entry/2020/09/20/205301
https://www.r-bloggers.com/2021/05/a-comparison-of-terra-and-stars-packages/
「自分が分析で使いたいパッケージに合わせて、ポリゴン・ラスタデータを扱うパッケージを選ぶと良い。」
tidyverseパッケージと相性がいいのはsfパッケージ & starsパッケージ (同じ作者です)
terraパッケージは従来のrasterパッケージの作者が開発したもので、ほぼ同じ機能が実装されているし、関数も豊富。
ややこしい作業をしないならばsf + starsでよいかも。
もし高度なラスタ演算が必要ならば、terraパッケージでラスタ演算–>starsで描画、という流れが良いかもね。
「自分が分析で使いたいパッケージに合わせて」ってどういうこと?
分析で使いたいパッケージの仕様書に書かれている関数のヘルプを読むのです。
例えば、景観分析でよく使うlandscapemetricsパッケージのcalculate_lsm()関数は、argumentsとして、Raster Layer, Stack, Brick, SpatRaster (terra), stars, or a list of rasterLayersを取れるので、terraでもstarsでも、rasterでも大丈夫。
(https://cran.r-project.org/web/packages/landscapemetrics/index.html)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
city_boundary <- sf::st_read(here::here('data',
'N03-20220101_27_GML',
'N03-22_27_220101.shp'))
## Reading layer `N03-22_27_220101' from data source
## `/Users/ch/Documents/GitHub/tm-data-viz/data/N03-20220101_27_GML/N03-22_27_220101.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 144 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 135.0913 ymin: 34.27182 xmax: 135.7466 ymax: 35.05129
## Geodetic CRS: JGD2011
class(city_boundary)
## [1] "sf" "data.frame"
summary(city_boundary)
## N03_001 N03_002 N03_003 N03_004
## Length:144 Length:144 Length:144 Length:144
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## N03_007 geometry
## Length:144 MULTIPOLYGON :144
## Class :character epsg:6668 : 0
## Mode :character +proj=long...: 0
glimpse(city_boundary)
## Rows: 144
## Columns: 6
## $ N03_001 <chr> "大阪府", "大阪府", "大阪府", "大阪府", "大阪府", "大阪府", "…
## $ N03_002 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ N03_003 <chr> NA, "大阪市", "大阪市", "大阪市", "大阪市", "大阪市", "大阪市…
## $ N03_004 <chr> NA, "都島区", "福島区", "此花区", "此花区", "此花区", "此花区…
## $ N03_007 <chr> NA, "27102", "27103", "27104", "27104", "27104", "27104", "27…
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((135.0918 34..., MULTIPOLYGON (((…
北摂の市町村のみ抽出してみる
hokusetsu <- c('豊中市', '池田市', '箕面市', '能勢町', '豊能町',
'吹田市', '高槻市', '茨木市', '摂津市', '島本町')
hokusetsu_city_boundary <- city_boundary %>%
filter(N03_004 %in% hokusetsu)
hokusetsu_city_boundary
## Simple feature collection with 13 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 135.3302 ymin: 34.7307 xmax: 135.6827 ymax: 35.05129
## Geodetic CRS: JGD2011
## First 10 features:
## N03_001 N03_002 N03_003 N03_004 N03_007 geometry
## 1 大阪府 <NA> <NA> 豊中市 27203 MULTIPOLYGON (((135.4478 34...
## 2 大阪府 <NA> <NA> 豊中市 27203 MULTIPOLYGON (((135.499 34....
## 3 大阪府 <NA> <NA> 池田市 27204 MULTIPOLYGON (((135.4392 34...
## 4 大阪府 <NA> <NA> 吹田市 27205 MULTIPOLYGON (((135.5068 34...
## 5 大阪府 <NA> <NA> 高槻市 27207 MULTIPOLYGON (((135.5838 34...
## 6 大阪府 <NA> <NA> 茨木市 27211 MULTIPOLYGON (((135.546 34....
## 7 大阪府 <NA> <NA> 茨木市 27211 MULTIPOLYGON (((135.5301 34...
## 8 大阪府 <NA> <NA> 箕面市 27220 MULTIPOLYGON (((135.4734 34...
## 9 大阪府 <NA> <NA> 摂津市 27224 MULTIPOLYGON (((135.5555 34...
## 10 大阪府 <NA> <NA> 摂津市 27224 MULTIPOLYGON (((135.5985 34...
ggplotで市町村名ごとに可視化してみる。日本語が含まれるので、https://github.com/Gedevan-Aleksizde/fontregisterer を参考にしてfontregistererパッケージをインストールしておく。
library(fontregisterer)
## You can see available font families by `names(quartzFonts())`
sans <- fontregisterer::get_standard_ja_fonts()['sans']
p_hokusetsu <- ggplot(data = hokusetsu_city_boundary,
mapping = aes(fill = N03_004)) +
geom_sf() +
scale_fill_grey() +
theme_minimal() +
labs(fill = '市町村') +
theme(text = element_text(family = sans),
axis.text.x = element_text(angle = 45, hjust = 1))
p_hokusetsu
まずはデータを読み込む。
library(stars)
## Loading required package: abind
lulc <- stars::read_stars(here::here('data', 'L03-b-14_5235', 'L03-b-14_5235.tif'))
lulc
## stars object with 2 dimensions and 1 attribute
## attribute(s):
## L03-b-14_5235.tif
## 50 :389732
## 70 :107739
## 10 : 57694
## 110 : 30196
## 100 : 18851
## 160 : 8083
## (Other): 27705
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 800 135 0.00125 JGD2000 FALSE NULL [x]
## y 1 800 35.3333 -0.000833333 JGD2000 FALSE NULL [y]
まずはプロットしてみる。
p_lulc <- ggplot() +
geom_stars(data = lulc) +
scale_fill_manual(values = c('yellow', 'beige', 'springgreen3', 'orange', 'red',
'grey30', 'grey70', 'brown', 'blue', 'lightyellow',
'lightblue', 'lightgreen', 'white'),
labels = c('田', 'その他の農用地', '森林', '荒地', '建物用地',
'道路', '鉄道', 'その他の用地', '河川地および湖沼',
'海浜', '海水域', 'ゴルフ場', '解析対象外')) +
theme_minimal() +
labs(fill = '土地利用') +
theme(text = element_text(family = sans),
axis.text.x = element_text(angle = 45, hjust = 1))
p_lulc
まず、市町村界を、rasterデータと同じJGD2000に投影変換する
hokusetsu_city_boundary_proj <- sf::st_transform(hokusetsu_city_boundary,
crs = sf::st_crs(lulc))
hokusetsu_city_boundary_proj
## Simple feature collection with 13 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 135.3302 ymin: 34.7307 xmax: 135.6827 ymax: 35.05129
## Geodetic CRS: JGD2000
## First 10 features:
## N03_001 N03_002 N03_003 N03_004 N03_007 geometry
## 1 大阪府 <NA> <NA> 豊中市 27203 MULTIPOLYGON (((135.4478 34...
## 2 大阪府 <NA> <NA> 豊中市 27203 MULTIPOLYGON (((135.499 34....
## 3 大阪府 <NA> <NA> 池田市 27204 MULTIPOLYGON (((135.4392 34...
## 4 大阪府 <NA> <NA> 吹田市 27205 MULTIPOLYGON (((135.5068 34...
## 5 大阪府 <NA> <NA> 高槻市 27207 MULTIPOLYGON (((135.5838 34...
## 6 大阪府 <NA> <NA> 茨木市 27211 MULTIPOLYGON (((135.546 34....
## 7 大阪府 <NA> <NA> 茨木市 27211 MULTIPOLYGON (((135.5301 34...
## 8 大阪府 <NA> <NA> 箕面市 27220 MULTIPOLYGON (((135.4734 34...
## 9 大阪府 <NA> <NA> 摂津市 27224 MULTIPOLYGON (((135.5555 34...
## 10 大阪府 <NA> <NA> 摂津市 27224 MULTIPOLYGON (((135.5985 34...
北摂の範囲だけの土地利用データにするために、lulcをclipする
hokusetsu_lulc <- sf::st_crop(lulc, hokusetsu_city_boundary_proj)
# ちょっと可視化したいだけならこれで十分
plot(hokusetsu_lulc)
では、オーバーレイしてみよう。
p_hokusetsu_lulc <- ggplot() +
# ラスタデータの描画
geom_stars(data = hokusetsu_lulc) +
scale_fill_manual(values = c('yellow', 'beige', 'springgreen3', 'orange', 'red',
'grey30', 'grey70', 'brown', 'blue', 'lightyellow',
'lightblue', 'lightgreen', 'white'),
labels = c('田', 'その他の農用地', '森林', '荒地', '建物用地',
'道路', '鉄道', 'その他の用地', '河川地および湖沼',
'海浜', '海水域', 'ゴルフ場', '解析対象外'),
na.value = NA) +
theme_minimal() +
# ベクタデータの描画
geom_sf(data = hokusetsu_city_boundary,
fill = NA, color = 'black', size = .5) +
# プロットの調整
theme_minimal() +
labs(fill = '土地利用') +
theme(text = element_text(family = sans),
axis.text.x = element_text(angle = 45, hjust = 1))
p_hokusetsu_lulc
## Warning: Removed 62938 rows containing missing values (geom_raster).
ラスタデータの各グリッドに、hokusetsu_city_boundary_projの市町村名の情報を空間結合する。
lulc_joined <- st_join(hokusetsu_lulc, hokusetsu_city_boundary_proj)
## Warning in st_join.stars(hokusetsu_lulc, hokusetsu_city_boundary_proj): st_join
## found 1508 1-to-n matches, taking the first match for each of those
空間結合した市町村名で各グリッドをグループ分けして、市町村ごとの土地利用を集計する。
lulc_joined_df <-
# ラスタデータ変換し、
# 1行に1グリッドの土地利用・市町村名が格納されたtibble形式のデータフレームを作成する
as.data.frame(lulc_joined) %>%
tibble() %>%
# 土地利用が格納された列名をlulcに変更
rename(lulc = L03.b.14_5235.tif) %>%
# 市町村名がNAではないグリッドだけ抽出
filter(!is.na(N03_004) & !is.na(lulc)) %>%
# 市町村ごとに土地利用クラスを集計
group_by(N03_004, lulc) %>%
summarise(N = n()) %>%
# 土地利用クラスのIDを名前に変換
mutate(lulcname = case_when(lulc == 10 ~ '田',
lulc == 20 ~ 'その他の農用地',
lulc == 50 ~ '森林',
lulc == 60 ~ '荒地',
lulc == 70 ~ '建物用地',
lulc == 91 ~ '道路',
lulc == 92 ~ '鉄道',
lulc == 100 ~ 'その他の用地',
lulc == 110 ~ '河川地および湖沼',
lulc == 140 ~ '海浜',
lulc == 150 ~ '海水域',
lulc == 160 ~ 'ゴルフ場',
lulc == 255 ~ '解析対象外',
TRUE ~ as.character(NA)))
## `summarise()` has grouped output by 'N03_004'. You can override using the
## `.groups` argument.
さあ、プロットしようか。
p_lulc_bar <- ggplot(data = lulc_joined_df,
mapping = aes(x = N, y = reorder(lulcname, N),
fill = lulc)) +
geom_bar(stat = 'identity') +
scale_fill_manual(values = c('yellow', 'beige', 'springgreen3', 'orange', 'red',
'grey30', 'grey70', 'brown', 'blue', 'lightyellow',
'lightblue', 'lightgreen', 'white'),
labels = c('田', 'その他の農用地', '森林', '荒地', '建物用地',
'道路', '鉄道', 'その他の用地', '河川地および湖沼',
'海浜', '海水域', 'ゴルフ場', '解析対象外'),
na.value = NA) +
theme_bw() +
facet_wrap(~N03_004) +
theme(text = element_text(family = sans)) +
labs(x = 'グリッド数', y = '土地利用区分') +
guides(fill = 'none')
p_lulc_bar
素直にgoogle検索 (“[package name] package” AND “[error message]”)
日本国内なら TokyoR に参加してみる https://tokyor.connpass.com/
@hagachi に聞いてみる
チュートリアル&あなた自身のデータでまずは手を動かしてみて、たくさんのエラーと最後の成功体験を楽しんでくださいな。